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Singular value decomposition for the 2D fan-beam 
Radon transform of tensor fields 
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Abstract — In this article we study the fan-beam Radon transform Dm. of symmet- 
rical solenoidal 2D tensor fields of arbitrary rank m in a unit disc D as the operator, 
acting from the object space L2(D;Sm) to the data space Z/2([0,27r) x [0,27r)). The 
orthogonal polynomial basis s^*^' of solenoidal tensor fields on the disc D was built 
with the help of Zernike polynomials and then a singular value decomposition (SVD) 
for the operator Cm was obtained. The inversion formula for the fan-beam tensor 
transform Drn follows from this decomposition. Thus obtained inversion formula can 
be used as a tomographic filter for splitting a known tensor field into potential and 
solenoidal parts. Numerical results are presented. 

1 Introduction 

The problem of determining vector or tensor field from the integral information 
arises in various applications, for instance in ultrasound probing of fluid or gas 
flows and deformed elastic media. In the flrst case it's required to determine 
the velocity vector field in the flow and in the second case — the stress tensor 
field. 

One of the most complete monograph of the tensor tomography is [23]. 
Reversibility and stability of different kinds of transforms of tensor fields on the 
Riemannian manifolds are studied there. In [5], [19] the solution of the vector 
tomography problem is reduced to the scalar Radon problem. An approximate 
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solution of the vector and tensor (of rank 2) tomography problem is given in 
[9], [10] with the help of polynomial non-orthogonal basis. 

More information and references about vector and tensor tomography prob- 
lems are given in [2], [22], [24], [25]. 

In this article we derive an inversion formula on the basis of singular value 
decomposition (SVD) for the fan-beam transform of tensor fields. To this end 
the orthonormal polynomial basis of solcnoidal tensor fields, supported in imit 
disk, are built from Zernike polynomials. In the scalar case thus obtained 
SVD corresponds to the known SVD for the Radon transform in the classical 
(parallel) formulation [6], [7], [12], [16], [15]. 

Unlike the scalar case, Radon transform of tensor fields has a non-zero kernel 
and it's possible to reconstruct uniquely (without additional information) only 
the solcnoidal part of a tensor field, so the inversion formula can be used as a 
tomographic filter for splitting a known tensor field into potential and solenoidal 
parts. 

This article is organized as follows: In Section 2 we formulate the problem 
of 2D tensor tomography. In Section 3 we review those part of the tensor fields 
theory that are needed in this paper. Section 4 contain a novel properties of 
Zernike polynomials. Sections 5 is devoted to the orthogonal polynomial basis 
in the space of solenoidal (divergence free) tensor fields and a singular value 
decomposition (SVD) for the tensor tomograph problem. A short description 
of the implementation issues and numerical tests are presented in Section 6. 

2 Formulation of the problem 

Let us consider the Cartesian coordinate system {x^,x'^) on the plane and 
let Tm denote for m = 0, 1, ... the space of all real- valued m-covariant tensors 

a := ai^...i^dx^^ ig) dx^^ ig) ... ig) da;*"* or a = {ai^...i^, h, ...,im = 1, 2}, 

where (g) is the tensor product and ftij .. ^^ are the components of a in the Carte- 
sian basis {x^,x^). Here and throughout we imply the summation convention. 
By Sm we denote the subspace of symmetric m-covariant tensors and there 
exists a canonical projection tr : — > (called symmetrization) onto this 
space defined by the equation 



where Ilm is the group of all permutations of degree m. A symmetric m- 
covariant tensor a = {ai^,,,i^,ii, ...,im = 1,2} has only m + 1 independent 
components which we denoted by afe, so that 




(2.1) 



«fc '■= "1...12...2' = 0, ...,m). 



(2.2) 
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Wc will always denote vector and tensor fields and any related quantities 
such as functional spaces by boldface characters. 

Let B := {{x^,x'^) € (^^1^2 (-^2^2 < 1} bg a unit disc on the plane M^. 

The symmetric m-covariant tensor field a(a;^,x^) defined on D can be treated 
as a mapping 

aiP^S™, a{x^,x'^) = {ai^...i^{x^,x'^), ii, = 1, 2}. 
The fan-beam Radon transform Vm of tensor field a(a;^, a;^) is defined by 

2cos(/3-<p) 

[Dma\{P,(p) := j e'^-e'^-...-e''^ai,...i^{cosP-lcos>fi,smP-lsmip)dl, 


(2.3) 

where /3 e [0,27r), 6 = (^') = , I/? - ^1 < ^ 



\ sin ip I 2 ' 

The difference between the parallel-beam and the fan-beam geometry is shown 
in figure 1. 



Figure 1: Left: parallel-beam scanning geometry. Middle: fan-beam scanning geometry. 
Right: an example of the fan-beam transform (the data function or sinogram) f(f3,ip), the 
angle fS defines the vertex point of the fan-beam projection /(/3, •) and the angle (p defines 
the direction of scanning. 



For 1/3 — </?| > we complete the definition of the fan-beam transform (2.3) 
with the condition 

[I?„a](/3,^) := {-ir+'[Vma]i(3,ip + 7T). (2.4) 
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Note, that the case m = corresponds to the fan-beam Radon transform 
Vq = 1? of a, scalar function a{x^,x^). 

Now, the problem is to recover the unknown tensor field a(a;^,a;^) in the 
unit disc D from the data function /(/?, (p), see Figure Ic, such that 

[V„^^i]{p,ip) = f{p,ip), {p,ip) e [0,27r) x [0,27r). 

This problem will be solved here by the SVD-method. 

3 Preliminaries 

In this section, we introduce the definition of SVD method and then review 
some facts from vector and tensor analysis [23] and, in particular, consider 
real- valued tensor fields in complex coordinates (variables) [26]. We define 
here some functional spaces of tensor fields — L2(D;Sto), for example, and 
also establish the notations that will be used in the sequel. 

3.1 Singuleir value decomposition (SVD) 

Now we define the concept of a singular value decomposition, sec [17], [18], [15], 
[21]. Let U and V be Hilbert spaces, and A be a compact linear operator from 
U to V, A £ £{U,V). Then there exists a sequence {(Tk}k>i of positive num- 
bers, monotonically tending to zero (or a finite sequence) and two orthonormal 
systems {uk}k>i C U, {vk}k>i C V, such that for all u G U we have a singular 
value decomposition 

oo 

Au = ^^ak{u,Uk)uVk, Auk = <JkVk, cri> a2> ... > 0. 

fc=i 

The adjoint of A is given by 

oo 

A*v = '^(Jk{v,Vk)vUk, A*Vk = (TkUk 
fe=l 

and the generalized inverse of A is 

oo 
fe=l 

Operator can be unbounded, so one can use a truncated SVD for its regu- 
larization 

/c<l/7 

where 7 is the parameter of regularization. SVD is one of the methods for solv- 
ing ill-posed problems and it allows to characterize the range of the operator, 
invert it and estimate an incorrectness of the corresponding inverse problem. 
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3.2 Tensor fields in complex coordinates 

Let's identify with the complex plane C by the usual way 



z := + ix^, = z:= x^ 



• 2 -2 

IX , 1 



-1. 



Let a = {aj^...i^ (x^, x^)} be an m-covariant real- valued tensor field in Cartesian 
coordinates (a;^, x^), then in complex coordinates or variables {z,z) it will have 
new components Ai^,,,i^{z,J), which are formally expressed by the covariant 
tensor law 

dx"^ dx^" 



-^il...im.{^^ ^) 

where the Jacobian matrix is 



and the inverse matrix of it is 



dz^^ " dz^^ 



(3.1) 



fdz'^ 


dz'^\ 


dx^ 


dx"^ 


dz^ 


dz^ 


\dx^ 


dx^J 



J-' = 



IS 




/dx^ 


dx^\ 


'dz^ 




dx^ 


dx"^ 


\dz^ 





1 1 

— i i 



d 


d 


1 / d 




d 


d 


1/9 . d 


dz^ 


dz 


2 \dx^ 


'dx^ 


dz"^ 


dz 


2 \dx^ ' 'aa;2 



Hciic; the formal partial derivatives with respect to 2:^ and z'^ are defined in 
the usual way 



(3.2) 



We shall write transformation (3.1) as 

a= {aii...i^(a;\a;^)} ^ A = {Ai^,„i^{z,z)}. 

From now on small letters will be used to denote tensor fields in the initial 

Cartesian coordinate system (a;^,.x^) and capital letters will be used for the 
same tensor fields in complex coordinates {z,T). 
An inverse relationship also takes place 



Oil.. .1^(2; ,X ) 



dz^'^ dz^ 



dx^^'dx^"^ * 
and we shall also write this as 

A = {Ai^...i^{z,z)} ^ a = {aii...i„(a;\a;^)}. 



(3.3) 



5 



A symmetric m-covariant tensor A could also be given by its components 



^fe ^1...12...2' {k = 0,...,m) 



(3.4) 



and subject to the conditions 

Ak=Am-k, {k = 0,...,m). 



(3.5) 



So we may image the symmetric tensor as pseudovector, expanding the one as 
a column array for convenience, that the following notations will be used 



/ am \ 

ai 



f Am \ 

Am- 1 

Ai 
V ^0 / 



(3.6) 



Taking into account the tensor law (3.1), (3.3) we get the formulae that hnk 
independent components (2.2) and (3.4) in pseudovectors (3.6) 



m — k k 

■m-k ™-k k 
r=0 ,s=0 

where A; = 0, 1, m and C] are binomial coefficients. 



(3.7) 
(3.8) 



3.3 Metric tensor G and pointwise inner product in com- 
plex coordinates 

On parity with covariant components of the tensor we shall also use its con- 
travariant components. In Cartesian coordinates {x^,x^) covariant and con- 
travariant components gij and g^^ of the metric tensor g are the same 



S={9ij} = {9''} 



1 
1 



(3.9) 



Thus contravariant components of the tensor a coincide with its correspond- 



ing covariant components, ai 



^n - Jm 'j'jje pointwise inner product (•,•) 



on Sm induced by the Euclidian metric g (3.9) is defined by the formula 

(a, b) :=a*i'-^™6,,i,...,„. 
In complex coordinates {z,'z) the metric tensor G has the following covariant 



{Gij} - I I 
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and contravariant components 

Contravariant components of tensor A in complex coordinates are obtained by 
raising indexes with contravariant components of the metric tensor (3.10) 

and the pointwise inner product of tensor fields is evaluated by formula 

/A — i»l»2---jm D. . . _ 4 . . f>ili2---im 

If tensors A — {Ak} and B = {B^} are considered as pseudovectors in 
complex coordinates then their pointwise inner product will be equal to 

m 

(A, B) = 2™ ^ Cl;^AkBm-k. (3.11) 

fe=0 

The pointwise norm of tensor A then will be 

m 

|A|2 = 2™^C^|Afe|2. (3.12) 

fe=0 

It is clear that the pointwise inner product is invariant, i.e. if a ^ A and 
b ^ B, than 

(a, b) = (A, B). (3.13) 
3.4 The space of integrable tensor fields L2(1D'; S^) 

Let L2(D;Sm) denote a Hilbc^rt space c;omprising real- valued symmetric m- 
covariant tensor fields on ID with the inner product, denoted by ((•, •)) 

((a, b)) = ((a, b))L,(„.s„) := // {a{x\x% h{x\x^))dV\ dV^ = dx^Adx^ 



and the finite norm 1 1 • 



I l|2 — II ||2 
l^ll — ll^llL2(D;Sm) 



:= jj{a{x^,x'^), a(a;\a;2))dy2. 

D 

In complex coordinates for a ^ A and b ^ B we have 
((A, B))= JJ{A{z,z), B{z,z))dV^ dV 



2 dz A dz 



-2i 

By virtue of invariance of inner product (3.13) the following equalities take 
place 

((a, b)) = ((A, B)), ||a|| = ||A||. 
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3.5 Differential operations on symmetric tensor fields 



We shall denote the class of real-valued m-covariant symmetric tensor fields 
a = {ai^___i^{x^ ,x^)}, whose all components are functions from (7*^(0), 1 < 
fc < 00 by C'^(D; S^). A subset of C'°(B; S„) whose finite support is contained 
in D will be denoted by Cg(D; S„). 

The operator of covariant differentiation V (in the vectorial case = grad) 

V:C°°(]D);S™)^C°°(D;S„,+i) 

in Cartesian coordinate system (x^^x^) is defined by equation 



Va := {ai^...i^-j} 



1,2 



The covariant differentiation V operates on any tensor field of rank m > and 
produces a tensor field that is one rank higher. For example, the gradient of a 
(co)vector field is a second rank tensor field. 
In complex coordinates we have 



VA = {Ai„„i^.j} 



J 



1,2 



where formal partial derivatives with respect to and are defined by (3.2). 
The operator of divergence 5 (in the vectorial case = div) 

5:C~(e;S™)^C=°(P;S„_i) 

in Cartesian coordinate system (x^,a;^) is defined by 



(5a := {a 



X3 



The divergence 6 can operate on any tensor field of rank m > 1 and above 
produces a tensor that is one rank lower. For example, the divergence of a 
second rank tensor field is a (co)vector field. 

In complex variables the divergence is calculated with the help of contravari- 
ant components G^^ of the metric tensor (3.10) 



SA 



5A- 



dz 



i2 _|_ 2^^n»2 



dz 



1,2 



}.3. 



14) 



A smooth tensor field a G C''{B; S^) is called solenoidal if its divergence equals 
to zero. The condition for the tensor field A to be solenoidal can be expressed 
in complex coordinates in terms of its independent components A^, ■■■,A^ 

( {Ao), + {Ai)^ =0 / Am\ 



= 

{Ak), + (Ak+i)^ = 
{A^_i), + {A^)^ =0 



or 



d_ 
dz 



A 

A2 

V ^1 J 



d_ 

dz 



A 

Ai 



= 0. (3.15) 
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The next differential operation on the symmetric tensor fields is the sym- 
metric inner differentiation d 

d:C~(D;S™_i)^C~(D;S„), 

defined in the following way 

d := ctV, 

where a is the symmetrization operator (2.1). 

A tensor field a e C°°(]D); S^) is called a smooth potential field, if for some 
tensor field v G Ci^(lD); S„i_i) with boundary condition v|aD = we have 
a = dv and v is the potential. 

The symmetric inner differentiation d in complex variables is calculated in 
the following manner. If a = rfv and a ^ A, v ^ V then A = dV and 

m — k dVh k dVh-i , , s . ^ 

Ak = + m>l, fc = 0,l,...,m . 3.16 

m oz m oz 

3.6 Orthogonal decomposition of the space L2(D; S^) into 
the sum of solenoidal and potential peirts 

Operators d and —5 arc formally conjugate and for a boimdcd region G with a 
piecewise-smooth boundary the Gauss-Ostrogradsky formula takes place 

j I {{dv, a) + (v, Sa)) dV^ = J (vv, a) dV\ (3.17) 

G dG 

where a € Sm, v G Sm_i arc smooth tensor fields, and u = {h'1,1'2} G Si is 
a unit covector of outward normal to the boundary dG, and is the operator 
of symmetric multiplication with the covector 1/ 

iv • ^ S^Ti-i-i, 

which is defined by the equation 

In terms of the Gauss-Ostrogradsky formula (3.17) we can define that a tensor 
field a € L2(D; S^) is solenoidal if the following equation takes place 

{{dv, a)) = jj {dv, a) dV'^ = (3.18) 

D 

for all smooth tensor fields v{x^ ,x'^) £ Co°(IQ'; S^-i)- 

We denote by H(D; Sm, 6) the graph space of 6 over L2(]D>; S^), i-e. 

H(P; S„, S) := {u € Ui^; S„) 5u e L2(P; S„)}. 
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It is a Hilbert space under the graph norm 

((u,v))h(D;S^,5) := ((u,v)) + {{Sv,Sv)), ||u||^(n.s^_5) := ||u||2 + ||5u|p. 

Finally, we define subspace of solenoidal tensor fields (i.e. which satisfies the 
equation (3.18)) 

H(D; S„, (5 = 0) := {a e H(P; S„, 5) | (5a = 0} 

and it is clear that this subspace is a completion of the set of smooth solenoidal 
tensor fields with respect to the norm || • || of L2(B;Sm). 

By Hjv(B;Sm,'5 = 0) we denote the finite-dimensional subspace of poly- 
nomial (of degree at most N) solenoidal m-covariant tensors fields. Then we 
have 

Ho(P; S„, S = 0)C Hi(©;„ , ^ = 0) C ... C H;v(B; S„, ^ = 0) C ... C L^i^; Sm[ 
and 

H(©;S„,^ = 0)=clos( [J HAr(B;S„,^ = 0)J , 

\N=0 J 

where clos means the closure in L2(D; S^). 

It is well known, see [8], [27], that a vector field can be represented as 
a sum of solenoidal and potential vector fields. The classical result in this 
direction belongs to H. Weyl and is connected with the decomposition of the 
L-i space of vector fields into the orthogonal sum of solenoidal and potential 
fields. The analogous result is true for tensor fields, see [11], [14], [23]. Namely, 
for u e L2(D; S^) we have 

u = a-|-rfv, ((a, dv = 0)), 

where a is a solenoidal tensor field and v € Hq(P; Sm-i)- Or, in another words 
the orthogonal decomposition 

L2(I]); S„) = H(D; S^, ,5 = 0) (iHi(D; S„_i) 

takes place, where the Sobolev space Hj(D; Sm_i) is a completion of the space 
of smooth tensor fields Co(D;Sto_i) with respect to the Sobolev norm || • jji, 
corresponding to the scalar product ((•, •))! that is defined by the formula 

((u, v))i = ((u, v)) + ((Vu, Vv)). 

3.7 Fan-beam Radon transform "D^ of tensor fields in 
complex variables 

Let's assume that some (constant) vector field is given in Cartesian coordinates 



9"^ ) ~ \sm.if 
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then according to the tensor law for contravariant coniponents its representa- 
tion in complex coordinates will look like 

Then we denote by theta"^ the tensor product 

iheta™ := O O ... O 6> = {6^' ■ 6^^ ■ ... • ^^-j, 

V V ' 

m 

and 0™ will be an m-contravariant tensor in Cartesian coordinates and in 
complex coordinates we have the tensor product 

0"* := ® ® ... (g) © = {6^1 • 0^=^ • ... • e^-}. 

^ V ' 

m 

It is clear that 0"* ^ 0*". As soon as the inner product of tensors is invariant 
(3.13), we get 

(a, 6>™) = aj^...jj^^ ■ 9^^ ■ ... • (9^'" = (A, 0™) = Aj,...j^Q^' ■ 0^' ■ ... • O^-". 

Thus we can evaluate the fan-beam transform (2.3) through the components 
of the tensor A{z,z) 

2cos(/3-¥') 

['Dma]{(3,ip) = J (0'",a(cos/3 — Zcos<^, sin/3 — Zsintp)) dZ 



t t 
= J (0", A(C,C))|dC|= J e^^..e^-A,,...,„(C,C)|dC|, (3.19) 

T{t,(p) T{t,(p) 

where t = e'^, T{t, (p) = -te^'^ ^eK-^v-fi)^ p g [0, 27r), (f G 
Here and in the sequel, we use notation 



...|dC| 



for a line integral along the line segment with end points zi, Z2 G K). 

At last we can get the fan-beam transform (3.19) in terms of components 



Ak and for ip e 



we have 



11 



[P„a](/?,^)= / ^C^e"^^e-(™-'=)'^A,(C,C)|dC| (3.20) 

■J f, n 



m „ m 

= 51^™^'^"'"'"''" / ^fc(C,C)|dC| =^C^e'(2'=-™)'^[DAfc](3.21 



= J2 (^^e'(— (3.22) 

fc=0 

Recall that for \P — (p\ > J the fan-beam transform VmSi is defined by 
condition (2.4). 

Now wc verify that the potential part of a tensor field is "invisible" for 
tensor transform T>m- Let a = c?v and a ^ A = dV, v ^ V. Substituting the 
potential tensor (3.16) in the expansion (3.20) and making evident evaluations, 
we get 



t ^ 



d\z\ 



„ /m—1 o-rr 



d\z\ 

k-1 dVk- 1 
dz 



d\z\ 



r(*,v) 

dVk ^ 



"/Afro -'^^^h ^--a.J'^i^i 

m—1 * 
m—1 

= E e'(''=-'"+'^^C^_i (y,(tj) - Vfe(r,r)) = 0. 



/c=0 

Here 

d i,„ d d 



f,f^- - - a- (3.23) 
a© oz oz 



is the derivative in the direction = ( ^2 1 = ( ^-iip 1 > written in the complex 
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form and after integration we take into account that the potential V vanishes 
on the boundary of the disc D. 

At the end of this section we resume, that we consider the operator "D^ as 
follows 

Vm : L2(©; S„0) ^ i2([0, 27r) x [0, 27r), 

and kerVm coincides with the space of potential fields dHo(D; 8^-1), so one 
can say that potential fields are "invisible" for the tensorial Radon transform 



4 Zernike polynomials 

We will identify complex plane C with as above. Let D = {z : |z|<l}be 
the open unit disc in C and ^2(10) denote a Hilbert space comprising square 
integrable (complex-valued) functions on D with the inner product denoted by 

((•,•)) 

{{a,b)) := JJ a{x\x'^)b{x\x-^)dV^, dV^ = dx^ A dx^ 



and the finite norm 



|a||2 :=JJ \a{x\x^)\^dV\ 



It is known that Zernike polynomials [4] form a complete orthogonal system 
(basis) over the Hilbert space iy2(D). Let z ^ reJ^ € 1$, r = \z\, tp = arg(z). 
Traditionally, Zernike polynomials, see [4], [16], [20], are defined by 

y„,i(r, V) := e"'^i?„,|;|(r), {-n < I < n, ), (4.1) 

where 

(n-m)/2 ^ 

D (j,) ._ V (-l)P — — r"-2P 

are the so-called real- valued Zernike radial polynomials [20] , defined for integers 
n and m so that < m < n and n — m is even. The family Rn,m{i") is related 
to the Jacobi polynomials 

Rn,m{r)=r"^P^^\2r'-l), 

2 

where the Jacobi polynomials [13] are given through the Rodriguez formula 
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for a,b > —1; n = 0, 1, 2, and s G R form a complete orthogonal system in 
the Hilbert space L2[—l, 1] of square integrable functions on [—1,1]. 

In this paper we will use another numbering of Zernike polynomials (4.1) 
and treat them as polynomials in z, z. For this reason the new notation Z^''' 
is introduced according to the transformation of indexes Z = n — 2fc in (4.1). 
Thus we get 

K,„-2fe(r, V) = ^"-^'^Pf l"-"=l^(2|z|2 - 1), (A: = 0, l,...,n). 
So, we define 



Z"''=(z,z) := (-l)'^K,„-2fc(r,V) = (-1)"^ 



k^n-2k p(0,|n-2fe|) 



(2|; 



1), (4.2) 



where n = 0, 1, 2, ... and A: = 0, 1, n. The first index (superscript) n indicates 
the degree of a polynomial Z"'*^ and the second superscript k denotes its order 
in a bunch \ Z"'". The multiplier (-1)'= in (4.2) was introduced 

for convenience of further computations. 
This definition can also be rewritten as 



E:=o C^fcC'^-fe^""'"^(l - zzYi-z)"-^ for k = 0,l, ... 



Z"'''(z,z) = 



(-1)"Z"'" '"{z,z) 



for k = 



where [•] denotes the integer part of a number. 
After the evaluation of (4.3) we get 



+ l,...,n, 
(4.3) 



2'"''=(z,z) 



E*^ ( T^k—sr^s r^k—s—k—s^n—k- 



for fc = 0, 1, 



n — k . 



for k 



+ 1, ...,n. 

(4.4) 

For example, the first Zernike polynomials up to the degree (order) n = A are 
zo.o = 1 

^1,0 = ^ ^1,1 = _j 

Z2.0 = ^2 ^2,1 = 1 _ Z2.2 = ^2 

2-3,0 = ^3 ^3,1 = 22 - 32% ^ 32^2 _ 22 = -z3 

^4,0 ^ ^4 24.1 ^ 3^2 _ 4^3^ ^4,2 ^ ^ _ g^^ + 622^2 Z^.S = Sz"^ - izz^ ^4,4 ^ ^4 



The Zernike polynomials are orthogonal in the unit disc D, obey the follow- 
ing orthogonality relation 



n + 1 



(4.5) 



14 



and their i2-iiornis are equal to 



n+ 1 



, (fc = 0, 



It allows the expansion of an arbitrary function a{z,z) G 1^2 (D) in terms of 
a unique combination of Zernike polynomials. 



(4.6) 



n=0 



fe=0 



Since we use the complex variables z and z and treat them as independent 
variables here, we'll sometimes write a{z) instead of a{z,z). The formal partial 
derivatives with respect to z and z are defined in the usual way by (3.2). 

In the following theorem we formulate in complex variables some novel 
properties of Zernike polynomials. 

Theorem 1. The following properties take place: 
(a) Zernike polynomials (4.3) have the differential representation 



1 a* 



z'~ \ z 



(n > 0, k = 0,1,. ..,n). (4.7) 



k\ dz'^ 

(b) Zernike polynomials (4.3) are the solution of the elliptic system 



(^"'"). = 

(^n,n)_ + (Z"'"-l)^ = 



+ (Z"''^-!)^ = 



(4.8) 



(Z".i)- + (Z"'% 



and satisfy boundary conditions 



Z"'^(t,i) = (-l)'=r-2'=, 1^1 = 1, (n>0, fc = 0,l,...,n). 



(4.9) 



(c) Zernike polynomials (4.3) can be represented in the form of Cauchy-type 
integral 



1 r t"(t-z)'' 

2^ J {t - z)''+^ 

ltl=l 



'Z"'''(z,z) forn> 0, k = 0,1,..., n 



dt= < 



(4.10) 



for n> 0, k > n or k < 0. 



Proof, (a) Let's first prove (4.7) for fc = 0, 1, 
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By Leibnitz formula 

^ . . (fc) 
z" \ z 



(fc) 



^n—k 



ELo Cfe«^'^v('=-'*) we get 

1 (fc) ^ r n (s) r 

{1-zz)'' (1 

-I 2 — ' L J z L 



s=0 



(k-s) 



J2 Ck{n - k){n -k- l)...(n -k-s + 1)^""'=- 



xfc(fc - l)...(s + 1)(1 - zzYi-zY 

k 

= klJ2CtC:,_,,z'^->'-^l-zzr{-zY-^ =k\Z^''^{z,z). (4.11) 



s=0 



Now let's substitute A; n — fc in (4.11). Taking into account that k < n — k, 
we get 



z" \ - — z 



i — k 



n—k 



s=0 
fe 



(1 - ZZ 



.-=\n—k 



z''{l - zzY'^ 

- {n—k—s) 



-1 (n-/c) 



= Y.'^^n-kk{k-l)...{k-s + l)z 



k—s 



x{n-k){n-k- l)...(s + 1)(1 - zzyi-z^-''- 
= (n - fc)! ^ C^C«_,(-l)"(^)"-'=-«(l - zzn-z) 



k—s 



-n.k 



(-l)"(n - k)\Z ' {z,z) = {n- ky.Z"'''-\z,z). 
And the assertion (a) follows. 

(b) It is easy to verify equations (4.8) by direct computation using formula (4.4). 
On the boundary of the unit disc D, due to the normalization p^°'l"~^'^l'(i) = i 
of Jacobi polynomials, we get 



Z"''=(iJ) = (-l)'=^"-2^ |i| = l, (n>0, k = 0,l,-,n). 
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(c) In (4.10) we take advantage of Newtonian binomial formula 



1*1=1 1*1=1 



di 



s=0 
fe 



1 /■ 

2^1 J {t - 0)'=+! 

1*1 = 1 



di 



^ '^^ ^ fc!27ri 



1*1=1 



(i-z) 



fe+i 



di 



1 



k — s 



1 d'' 



A:! dz'^ 

k' 



1 

k\ dz^ 



k 
,s=0 



Z I z 



Theorem 1 is proved. □ 



4.1 Fan-beam Radon transform of Zernike polynomials 

The fan-beam Radon transform "D of a scalar function a{x^,x'^) is defined by 
(3.19) for m = 0. We have 



[Va]{(3,^)= J a(C,C)|dC|, 



(4.12) 



where t = e'"^ , (i e [0, 27r), r(t, f ) -te^''^, |/3 - (^| < = arg(t - r). For 

\(3 — (p\ > ^ we complete the definition of the fan-beam transform (4.12) with 
the condition 

[Da]{l3,ip) = -[Va\{p,<p + n). 

Theorem 2. The fan-beam Radon transform VZ^''^ of Zernike polynomials 
Zn,k equaJs to 



[vz-'^m^) 



2gi{n-2k)ip 

n + 1 



cos[(n + l)(/3 — (^)] for n = even 



(4.13) 



isin[{n + 1){(3 — (f)] for n = odd, 
where (3 e [0, 27r), e [0, 27r). 

Proof. At first we introduce the auxiliary polynomials X"'*' defined by 



X"''=(2,^) := 



1 9*^-1 



, (n > 1, fc = 1, n). 
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Then the next equations follows directly from (4.7) 

= = -Z^^'^-\z,z). (4.14) 
Using (4.7) we can verify, that 

k 

Then combining above and (4.9) we obtain the boundary conditions 

X"''=(tJ) =0, |i| = 1. (4.15) 

Let's compute the fan-beam transformation of Zernike polynomials. To this 
end we use previously obtained derivatives (4.14) for computing by (3.23) the 

derivative of X"'*' in the direction = ( 2 1 = I -i(n ) . So we get 



/iY""')fc fiVf^yk fiV^^f^ 

a© oz oz 

or the same in another form 



Z'''''{z,z) = e-2i¥'Z".'=-i + e"'*'- 



This equation combined with (4.15) is used for evaluation of the next integral 
t t 
y"2-"''=|dC| = e-2'^y"z"''=^^VCI +e"''^(X"^'^(tJ)-X"''^(T,T)) 

T T 

t 

= J Z"''=-i|dC|, 



where t = e'^, T(t, (p) = —te^^"^. Unwrapping the recurrence relation gives 
t t 



j Z"''=|dC| = e-^*^'"^ j Z"'0|dCI- 



The last integral is computed directly, taking into account that Z^'^{z, z) = z^ 

t t \t-T\ 

J Z"'0|dC|= J C"|dCl= I (T + se'*')"> 





-(r + se'^)"+^e-'^ '* ^' = e-'^(r+^ - r"+^^ 

1^ ' n + 1 ^ 
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Finally, we get 



2gi(n-2fc)¥' 



n + 1 



n + 

cos[(n + l)(/3 — for n = even 
i sin[(n + !)(/? — y>)] for n = odd. 



Theorem 2 is proof. □ 



5 Construction of the orthogonal polynomial 
basis and SVD 

In this section we describe the construction of orthogonal polynomial basis in 
the space of solcnoidal (divergence free) tensor fields H(D; S^, 8 = 0). 

We are given a polynomial of degree N solenoidal m-covariant tensor field 
a e Hjv(D; Sm, (5 = 0) and in complex coordinates we have a ^ A = {Ak}. As 
was mentioned earlier, we use a pseudovector notation for the tensor field A 



A = 



The condition (3.5) now looks like 



( Ara\ 
Am — 

V ^0 j 



^ Am 


\ 


A 


1 


Ai 




\ Aq 


/ 



(5.1) 



/ ^0 \ 

Ai 

Am — l 
\Am J 



(5.2) 



For m, n > 0, fc = 0, n + m and 2fc 7^ to + n we define polynomial of degree 
n symmetric tensor fields (in complex variables ) 



s^r' ■■= (-1)" 



^n,k—m+l _|_ 



,,/c-l 



V 



7n,k—m 



o(-m) 



\ 



Zn,k-1 _^,k-m+l 



.,fc-l 



n.k—m 



-z" 



where for convenience we set Z"'*^ = for fc < or fc > n. 



(5.3) 
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For 2fc = m + n wc have two cases: 
the first one is when m and n = even 



j(+m) 



the second one is when m and n = odd 



, s(-r]„ := 0; 



0, s 



(-m) 



(5.4) 



(5.5) 



Remark 1. Polynomial tensor fields (5.4) and (5.5) can be evaluated by 
general formulae (5.3), but then the result should be divided by 2. 

Remark 2. Note, that for + s = m + n the equation 

s?r^ = (-i)"s(±™) 

takes place, therefore in (5.3) we set only fc = 0, 1, [^^^^j^] , where [•] defines 
the integer part of a number. 

So, if wc now make transformations (3.3) or (3.7) from complex variables 

to real variables 



,(±m) 



n + m 



(5.6) 



then we get polynomial real- valued tensors s'"^^^ in Cartesian variables {x^,x'^). 

Lemma 1. The tensor Gelds s^^"^ (n = 0, TV, /c = 0, 1, [^^] ) 
defined by (5.6) form an ortliogonal basis of finite-dimensional subspace 
Hjv(P;S™,5 = 0), thus 



dimHjv(D;S„,5 = 0) = 



(iV+l)(iV- 



2m) 



Proof. Consider a tensor field a G H7v(lD'; S^, (5 = 0) and let a ^ A. 
Expand each component A/, of pseudovector (5.1) in the sum of Zernike poly- 
nomials and use the solenoidality condition (3.15). Taking into account the 
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property (b) from Theorem 1, we get an expansion of tensor A into the sum 



■^m—l 

\ Ao J 



N n-\-m 

= E E 
71=0 fc=0 



^n,k—m-\-l 



Prom another hand, (5.2) yields 



Ai 

Am — l 

\A^ J 



N n+m 

Z] c„,fc 

n=0 fe=0 



— n,fe— m+1 

Zj 



(5.7) 



n=0 /c=0 



— n,n — fe+1 

Zj 



-=7i,n— fc+m— 1 
Z 

-=n,n—k-\-m 



^n,fe— 1 



m+l 

Comparing the last expression with (5.7), we obtain 

Cn,k = {-'i-TCn,n+m-k, (fc = 0, n + m). (5.9) 

Splitting the coefficients Cn,k in (5.7) into the real and imaginary parts 
_ -u I I. — n fn + m 



taking into account (5.9) and definition of S^™\ we finally get 

in,k-\ _^-^ri,k-m+l 



( A^\ 

Am-l 



Ai 
V ^0 ) 



N 



n=0 fe=0 



7n,k—m 



( Z"^^^ - z" 

2n,k-l _-^n,fe-TO+l 



+ '^K,k 



2^n,k-m+l _ ■^"''^ ^ 
y Zn,k-m _ J 



N 



= E E {-iran,kS^r:r\^,^)-bn,kSij\z,z). 



n=0 fc=0 



The sign * here means that in the case of even n and even m, the coefficient 
b„ n+m should be set to 0, and in the case of odd n and odd m the coefficient 

"> 2 

a„ n+m should be set to 0. 
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So, we have that the tensor field a is a linear combination of polynomial 
tensor fields (5.6). 

Now we show that polynomial tensor fields (5.6) are orthogonal. Let k ^ 
s, k, s = 0, 1, [^^^^^] and remark, that then k + s ^ m + n take place. Using 
formula (3.11) wc have 

= ((s?r\ s(±-))) 

Taking into account the orthogonality of Zernike polynomials we obtain 

s(±'")))=0, k^s. 

Let's now evaluate the norms of polynomial tensor s^^^'' by the formula (3.12). 
At first we consider the case n + m ^ 2fc, then 



s. 



(±m)||2 
'n,fe II 



P P m 



■n,fe-m+p _|_ ■^.fe-P|2 



p=0 
m 

p=0 

fe 

e^m+l ^ ^ ^p y^^n,k—p 

p=k—n 

(m) 



,k—p I 



om+l 



(5.10) 



where coefficients are defined by 



p—k—n 



(5.11) 



Remark 3. In this formula for convenience we set = if p < or 

p > m. 

If n + m = 2k, i.e. k = then taking into account the above calculations, 
we get 

for n=odd 



,(+™) ||2 



2^77 

. n + l' 



Am) 



for n=even, 
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|„(-"*) ||2 
P„ "'+"11 



= < 



Obviously, 



2"7r 

( 

n + 1 



m+N 



„(m) 



for n=odd 



for n=even. 



Lemma 1 is proved. □ 

Lemma 1 and the definition of the subspace of solenoidal tensor fields 
H(e; 8^,(5 = 0) yield 

Corollary 1. Polynomial tensor fields s^^™^ (n > 0, A; = 0, 1, [^] ) 
form an orthogonal basis in the subspace of solenoidal tensor fields 
H(D;S„,^ = 0) CL2(D;S™). 

Now we are in a position to define the SVD for the fan-beam Radon trans- 
form Dm- 

Theorem 3. The singular values of the operator (2.3) 
Vm : L2(D;S„) ^ L2([0,27r) x [0,27r)) 

are given by 



(m) (±m) 



(n+l)2 



n + m 



where coefficients a^™-* are defined by the formula (5.11). If a solenoidal real- 
valued symmetrical tensor field a(a;^,a:^) e L2(]D>; S^) has an expansion 



a{x\x'^) = X! Z] (±m)|| («».fc4t^'(^^^^) + V/=si7r^(a;\a;2)) , 

n=0 fc=0 IPn,fc II 

(5.12) 

then the fan-beam Radon transform PrnS has the following singular value de- 
composition 

rn+mi 
OO I — 2 — i * 

M(/3,<^) = 5^ 4? («n,fe/^rH/5''^)+^".'^/nrH/3,<^)), (5-13) 

n=0 k=() 

where singular functions arc 

'cos[(n + !)(/? — ip)] cos[(n — 2k + m,)ip] 



_sin[(n + l)(/3 — </?)] sin[(n — 2k + m)ip], 
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cos[(n + !)(/? — if)] sin[(n — 2fc + m)(p] 



^sin[(n + l)(/3 — (p)] cos[(n — 2k + m)(p], 
when n > 0, A; = 0, [^^] and2k^m + n and 

1 j cos[{n + - ip)] 



V2n I 0, 







\/27r lsin[(n + l)(/3-(p)], 



when 2k — m + n. In all expressions above top line corresponds to the even 
values of n, and bottom line — to the odd n. The sign * in (5.12) and (5.13) 
near by the inner sum denotes that in the case of even n and m the coefhcient 
b„ n+m should be set to 0, and in the case of odd n and m — the coefficient 
g^^ n+m respectively. 

Proof. Note, that the system of functions //^^J forn > 0, A: = 0, 1, [^^±21] 
is the subsystem of the standard orthonormal basis of i2([0, 2tt) x [0, 2tt)) and 
there is the basis of the image of the tensor fan beam transform Vm- 

Let's evaluate now the fan-beam transform Vm for basis polynomial tensor 
(5.6). For this we introduce 



A := 



m+1 
^n.k—m 



I 



A„_, = Z"'*^-" and B := 



Hence in case of n + m 7^ 2A; we get 



-p^,k—m-\-l 
Zj 



-y=n ,k — s 
■ Zj 



si"!r^ = (-ir(A + B), siT^ = i(A-B). 
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Then by using formulae (3.21) and Theorem 2 we have 



2gi(m-»-2;/)v; 



m+1 



s=0 



n + 



5=fc— n 



n + 1 



cos[(n + — <^)] for n = even 

i sin[(n + l)(/3 — if)] for n = odd 
cos[(n + — <^)] for n = even 

^ i sin[(n + l)(/3 — <^)] for n = odd, 



where are defined by the formula (5.11). 

Analogically, using the formula (3.22) and Theorem 2 we get 



[VmB] = V„ 



-=n,fe— m+1 



— /i , k 

\ z' J 



s=0 



Op — i(m+n— 2fc)(^ ^ 



n+ 1 



X < 



s=k—n 



cos[(n + l)(/3 — (^)] for n = even 



_ — i sin[(n + l)(/3 — <^)] for n = odd 



2g-i(m+n-2fe)v3 

n+1 



cos[(n + 1)(/? — <^)] for n = even 
— isin[(n + l)(/3 — <^)] for n = odd. 
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Prom two formulas, derived above, it follows that 



{0M = 



(+m) 



2(-ir 

n+ 1 



4a, 



(to) 

n,k 



(/3,<^) = (A + B)] 

'cos[(n+l)(/3-^)] 

.Q,(™)('gi("-2fe+m)¥> _|_ g-i(n-2fe+m)v3^ X < 

[isin[(n + l)(/3-(p)] 
cos[(n — 2k + m)(^] cos[(n + — (p)] 



n+ 1 



X < 



^sin[(n — 2k + m)(f] sin[(n + — ip)] 



(+m)ll2 I cos[(n — 2fc + m)(/?] cos[(n + l)(/3 — (p)] for n=even 



sin[(n — 2A; + m)<p] sin[(n + !)(/? — 1^)] for n=odd. 



By the same way we evaluate the fan-beam transform of the other part of basis 
for n + TO 7^ 2k. 



™ n.k 



™ n.k 



(/3,<^) = i[D„(A-B)] 



2a, 



(m) 



^/ 2fc+m)¥' -i- — i(n— 2A;+m 



n + 



■cos[(n + l)(/3-(p)] 



4a, 



[sin[(n+l)(/3-^)] 
(m) [ sin[(n — 2fc + to)<^] cos[(n + l)(/3 — <^)] 



n.k 



n+l 



(-m)||2 



cos[(n — 2k + rnjLp] sin[(n + !)(/? — 

sin[(n — 2k + m)ip\ cos[(n + l)(/3 — tp)] for n=even 



7r2'"+i 



X < 



cos[(n — 2fc + to)^?] sin[(n + l)(/3 — (^)] for n=odd. 



Consider the case n + to = 2A; and n, m = even, then 



4iisl+r]ji^ 

2 

7r2'"+i 



(+m) 



4a^"l^ 

fa \ n.k 



cos[(n+l)(/3-^)] 



2(n + l) 



cos[(n + — ip)] for n=even 







for n=odd. 



If n + m = 2fc and n, m = odd, then 



(/3,<P) = 



■'^m'=>„ m+n 



(/3,V) = 



4a, 



(m) 
n,k 



2{n+l) \sm[{n + l){P-(p)] 
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4||s 



( - ) 112 



for n=even 

sin[(n + !)(/? — y)] for n=odd. 



Using equations for norms (5.10), we get (5.13). Theorem 3 is proved. □ 
At the end of this section we present some examples. 

Example 1. Let's take for instance m = 0, that corresponds to the scalar 
field, hence we have a{x^, x'^) = A{z,z) and the orthogonal basis s'"^^^ = S^^^^ 
in L2(ID)) is 



.(+0) _ 



'(-l)"2Re Z 



= < 



n.k 



„(+0) _ yn,^ 




for 2k 

for 2k = n, 

for 2k ^ n 

for 2k = n. 



where n > 0, k = 0,1,..., [f] . We have dimHjv(©; So, (5 = 0) = (^+iK^+2) 
and singular values are 



Stt 



p (n>0, fc = 0,..., [|]) 



_ (±0) 



Example 2. For m = 1 one gets covector field sl{x^,x'^) = {01,02}, which 
in the complex variables according to the tensor law has the representation 

A(z, z) = I ^ , ^ I • Dimension of the finite-dimensional subspace 
HAr(D; S^, 5 = 0) equals to ^nd singular values are 



J±i) 

'n,k 



Alt 

n+1 
Stt 

n + 1 



if n > and k = 



if n > 1 and k = 1, .... 



n+1 



The polynomial orthogonal basis s^"*"^^^ of the space of solenoidal covectors fields 
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n,k n,k 



n,k n,k 



H(©; Si, ^ = 0) looks as follows 

' (-1)" {Z"''^ + T"''~\ Z"''^-! + Z"'*"} for 2A; ^ n + 1 

S^+^ = {0, 0} for 2k = n + l, 

1 l^n.fc _ Z"''^-! - Z"''} for 2fc 9^ n + 1 

= Y {^"'"^'^"'"^"^l for2fc = n+l, 

where n > 0, fc = 0, 1, [^] . 

Exeimple 3. For m = 2 we have a symmetric second-order 2D tensor field 



A{z,z) 



an 0,12 

021 £122 

An Ai2 
A21 A22 



, which in complex variables has components 
ail - 022 - 2iai2 an + 0122 



4 

an + 022 



ail - £122 + 2iai2 



We also have inverse equalities 
Singular values for this case are 



2(^12 + Re All) -2ImAii 
-2ImAii 2(^12 -Re All) 



J±2) 



2tt 
n + 1 

n + 1 

n + 1 

Stt 



if n > and k = Q 



if n = and k = 1 



if n > 1 and k = 1 



n+21 



- if n > 2 and A; = 2,..., 
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where n > 0, fc = 0, 1, [^^] and basis tensor fields are 



,(+2) 
'n.fe 



j(+2) _ 



n,k 



= < 




^n,k-l 



j(-2) _ 



n.k 



= < 



^n,k-2 



1 ( _ 

for 2fc 7^ n + 2 



c(-2) 

"-5 71+2 




:n,k—l 



7n,k—2 



^n^k 



Also in this case we have dimHjv(ID'; S2, (5 = 0) = _ 



6 Implementation 

Scalar and vector cases of the inversion formula were numerically implemented 
and tested. The algorithm consists of 3 parts: solving the direct problem (that 
emulates the data acquisition in real life), finding coefiicients of the polynomial 
that represents a function being reconstructed and evaluation of this polynomial 
on a grid for visualization. 

In the scalar case, given a test function, defined by its values on a rectan- 
gular grid, the direct problem was solved by computing integrals (4.12) for the 
number of discrete values dp. Lpq 



[Da]{l3p,pp--+ipg)=fp,g, {p,q = 0,1,.. .,M+1). 



(6.1) 



Bilinear interpolation was used to get the values of the original function between 
knots. So the obtained data set is an (M + 2) x (M + 2) matrix of (/p,g) values 
that serves as an input for the inversion algorithm. Consider the scalar case 
for instance. Then, the function a{x, y) is approximated by the polynomial of 
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degree N < M (note, that in this section we use notations x = and y = x'^) 

N [ri/2]* 

aN{x,y) = 2^ ^ [^(a„,fcCos((n-2fc)^)-6„,fcsin((n-2A;)V')) 

n=0 fe=0 

X(_l)fe^n-2fcp^(0,n-2fc)^2r2 - 1)1 . (6.2) 

Here (r, ^) are the polar coordinates of the point (a;, y) and the sign * means 
that in the case of even n the coefficient a„ [„/2] should be divided by 2 and 
^n,[n/2] should be set to 0. Then the fan-beam transform of (6.2) will look like 

^ 4 '"Z^* /cos[(n+l)(/?-^)]cos(n-2fc)(p 
^ [sm[(n + l)(<^-/?)]sm(n-2fc)(^ 

X + ^^('^ - - 2^)^ (6.3) 

lsin[(n + !)(/? — (/?)] cos(n — 2A:)(/5, 

where the upper lines in the braces are used for the even n and the lower lines 
— for the odd n. The sign * means the same as in (6.2). After the substitution 
of (6.3) into the (6.1) we get a system of linear equations for determining 

(jV + -^-H-^ ^) unknown coefficients a„,fe and hn,h- 

2-K 



In the case of regular scanning scheme j3p = pe, ipg = q-, e 



'2' M + 2 

an explicit formulas for determining coefficients an,k and bn,k were derived, 
provided that M = N 

^ M+lM+l 

a„,fe = (-l)'- Y,fp,<,^^[e{pi2k-n) + ^{2k + l))], (6.4) 



M+1 M+1 
^ p=0 q=0 

Analogical formulae for the parallel-beam geometry can be found in [16]. The 
implementation of formulas (6.4) and (6.5) uses FFT and requires 0{N'^ log2 A'') 
operations, sec [3]. 

After the coefhcients a„,fc and bn,k are found, the polynomial that represents 
the reconstructed function is effectively evaluated using a recurrent formula, 
see [3]. 

The inversion algorithm was also tested under the presence of a noise in the 
input data (sinogram). Uniform and Poisson random distributions were used 
for this purpose. 

A representative set of numerical tests was performed for scalar and vector 
cases of the inversion algorithm. Some of the results are shown on figures 2-4. 
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On the top of the figure 2 there are original function (to be reconstructed) 
on the left-hand side and its sinogram (an input data set for the inversion 
algorithm) on the right-hand side. The middle row contains reconstructions 
from 32 and 256 fan-beam projections (free of noise). The number of terms 
in SVD was 30 and 254 respectively. The bottom row contains examples of 
reconstruction from noisy data. A random noise was superimposed on the 
sinogram. The _L2-norm of the noise was 10% of the L2-norm of the sinogram. 
1024 noisy fan-beam projections were used. The number of terms in SVD that 
were taken for reconstruction are 1022 and 254 respectively. It's possible to 
reduce the noise in the output image by taking less terms in the SVD. 

Another example of scalar tomography is shown on the figure 3. Again, 
the original unknown function (the fast oscillating one, with fine features) is at 
the top row, on the left-hand side and it's sinogram is on the right-hand side. 
The middle row contains reconstructions from 32 and 512 fan-beam projections 
(free of noise). The number of terms in SVD was 30 and 510 respectively. The 
bottom row contains examples of reconstruction from noisy data. A random 
noise was superimposed on the sinogram. TheL2-iiorm of the noise was 10% of 
the L2-norm of the sinogram. 2048 noisy fan-beam projections were used. The 
number of terms in SVD that were taken for reconstruction are 1022 and 510 
respectively. Again, one can observe significant enhancement of reconstruction 
when only part of terms are taken in SVD. 

The figure 4 illustrates the vector case of the inversion algorithm. The first 
solenoidal vector field (the top row, where the first component ai is on the left 
and the second component 02 is on the right) is defined by the formulae 

ai{x, y) = 2xycos{x^ + y^) + cos(6a;y) — &xys,m{Qxy), 

a2{x, y) = — sin(a;^ + y^) — 2x^ cos{x^ + y^) + 6y^ sin(6a;y). (6.6) 

Another vector field (the middle row) was obtained from the previous solenoidal 
vector field by adding the potential part 

ai {x, y) = 2xy cos(.x^ + J/^) + cos(6a;y) — 6xy sm{6xy) 

+ 2'KX cos(7r(.T^ + y^)), 

a2{x, y) = — sin(.x^ + y^) — 2x^ cos{x^ -\- y^) + 6y^ sm{6xy) 

27rycos(7r(x'2 + y^)). (6.7) 

As it can be seen (on the bottom row of the figure), reconstruction from these 
two vector fields is identical and contains only the solenoidal part of the vector 
field. 

The figure 5 illustrates another solenoidal vector field (the top row). Its 

reconstruction from the 20 irregular fan-projections is in the middle row. Here 
the positions of the fan-projection centers are shown as white dots on the 
boundary of the circle. Scanning was performed only over those lines, whose 
endpoints belong to this set of 20 points. The last reconstruction (the bottom 
row) was made under the presence of noise in the sinogram. The noise level 
was 3% (again, in L2-norm). 
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Figure 2: Top row: image of the original function (on the left) and its sinogram (on the 
right). Middle and bottom rows: reconstructions of the function from different number of 
fan-projections and under the presence of noise in the sinogram, see the text for detailed 
explanation. 
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Figure 3: Top row: image of the original function (on the left) and its sinogram (on the 
right). Middle and bottom rows: reconstructions of the function from different number of 
fan-projections and under the presence of noise in the sinogram, see the text for detailed 
explanation. 
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Figure 4: Top row: solenoidal vector field (first component on the left, second component on 
the right) given by formulae (6.6). Middle row: non-solenoidal vector field given by formulae 
(6.7). Bottom row: reconstruction from 20 fan-projections is identical in both cases, relative 
error 0.21%. 
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Figure 5: Top row: original solenoidal vector field. Middle row: its reconstruction from 

20 irregular fan-projections, relative error 2.6%. Bottom row: its reconstruction from the 20 
regular fan-projections under the presence of noise (3%) in the sinogram, relative error 1.6%. 
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7 Conclusion 



The novel inversion algorithm for the tensor tomography problem was devel- 
oped and mimcrically implemented. The algorithm is based on the SVD of the 
tensor Radon transform that allows to characterize the range of the operator, 
invert it and estimate an incorrectness of the corresponding inverse problem. 
The algorithm can also be used with noisy measurements. 
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